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In this paper there is presented method for ab initio calculation of the phonon spectra. The method 
is based upon a direct calculation of the dynamical matrix via second derivatives of the total energy. 
The pseudopotential technique in plane- wave basis set was used to calculate the total energy within 
the local density approximation (LDA) . For the change of the electron density there is derived the 
self-consistent equation which is solved analytically with no use of iterations. In this paper the 
attention is paid only to non-metallic systems. 
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I. INTRODUCTION 

Increasing accuracy of methods for calculation of the 
electronic structure of solids enable us to obtain the vari- 
ety of properties of the solid state. In this paper, there is 
described new scheme for ab initio calculations of phonon 
spectra in non-metallic systems. 

The problem of phonons as the perturbation of the per- 
fect crystal lattice is used to be studied in diffcrgiit ways. 
The first one is the so-called direct approachau based 
upon the comparison of the properties of unperturbed 
and perturbed states. The perturbed state is solved with 
use of the supercell geometry. This method enable us 
to study both the linear and nonlinear properties of the 
crystal lattice. Disadvantage of this method is that for 
each phonon wavevector q we need a corresponding su- 
percell. Only a few phonons can be described by the su- 
percell containing only a few cells, and this is the reason 
why the use of this method is limited only to few special 
long-wave phonons. Thus, this supercell technique seems 
to be insufficient as there are a lot of quantities, which 
must be computed by integrating over phonon frequen- 
cies of the whole Brillouin zone. 

The second way is to calculate the dynamical matrix 
which can describe dynamical properties of solids with 
use of the linear response methods. 

In this paper there is presented method for ab initio 
calculation of phonon spectra of non- metallic systems. At 
first the self-consistent electron charge density is calcu- 
lated. Then, the self-consistent variation of charge den- 
sity for an arbitrary phonon wave-length is calculated. 
After that, that one is used for calculating of the dy- 
namical matrix. The presented scheme is restricted to 
non-metallic systems because the variation of the num- 
bers of states with respect to perturbation in non- metallic 
systems is almost zero. The situation in metals is more 
complex. 

The presented method differs from the other Jinear 
response methods, they, dielectric matrix methodo, the 
method of Baroni et ala, the method of King-Smith and 
NeedsO and the method of Rignanese, Michenaud and 



GonzdZl'l. The presented method is based on direct calcu- 
lation of the dynamical matrix via the second derivatives 
of the total potential energy of the crystal lattice. 
Throughout the paper, the hartree units are used. 



II. DYNAMICAL MATRIX 

In the harmonic approximation the theory of lattice 
dynanucs can be formulated in terms of the dynamical 
matrixu. Then it can be written 



[D{q)^i,^q)]W{(l) = 0, 



(1) 



where D stands for the dynamical matrix, q is the 
wavevector of the phonon, to is the frequency of the 
phonon and W is the eigenvector of the dynamical ma- 
trix which represents the polarization of the phonon. The 
dimension of the dynamical matrix is 3N x 3N, where N 
stands for the number of atoms per ceU.. Elements of the 
dynamical matrix can be expressed asH 



Dai3{q,fiiy) = 



1 



v/Mpw: 



^ Aajsih, ixv) exp (-iq • h). 



(2) 



where M^, M^ are masses of atom /i and v^ a and [3 are 
Cartesian directions and h is vector of the translation 
symmetry. Elements of the matrix A are defined as 



Aap{h.,^iv) 



def 



B^Vl 



dLL,0 o f.h 



(3) 



where Vl stands for the total potential energy of the 
crystal lattice, u^'" is a vector giving displacement of 
atom ^ away from equilibrium in direction of a axes. The 
u o' is a vector of displacement away from equilibrium in 
the direction of (3 axes of atom v in elementary cell with 
the translation vector h from the cell containing atom /i. 
Defining the operator 









^ h 



OK''' 



u=0 



one can write the dynamical matrix as 



(4) 



(5) 



where asterisk means the complex conjugated operator. 
To solve the Eq. (|l|) we have to know elements of the 
dynamical matrix. They can be calculated on the basis 
of elements of matrix A. To calculate elements of the 
matrix A we need the second derivatives of energy Vl for 
which it can be written 



Vl = Vew + Eel , 



(6) 



where Vew is the Ewald ion-ion interaction energy and 
Eel stands for energy of the set of electrons. 

To find the contribution of the Vew to the dynamical 
matrix is easy because the analytical expressions for that 
one are well knownlld and therefore in the rest of this 
paper we study only term E^i . 



III. THE ELECTRON SYSTEM ENERGY 

Using the conventional density functional formalismEJ 
in the pseudopotential franaework, the total energy of the 
set of electrons is given bya 



Eel=E,-]^j p{v)Ve{v)d\ 



(7) 



where p is the charge density, Vc is the Coulombic po- 
tential, pxc is the exchange and correlation potential and 
Exc stands for the exchange-correlation energy. The E^ 
is the sum of eigenvalues 



Ee = y^^njEi, 



(8) 



where eii-are eigenvalues of the one-particle Kohn-Sham 
equation^ 



(T +U + Ve+ Pxc)A = £^'^^ 



(9) 



^i are self-consistent eigenfunctions, T is operator of ki- 
netic energy and U stands for pseudopotentials of the 
ionic sites. The exchange-correlation potential p^c is 
given in the local density approximation by 



Pxc{v) 



5p{v) 5p{r) 



(10) 



The variation of E^i can be then expressed as 

5Eei = <S^n,£, - ]^5 j p{Y)Ve{v)d\ 

i 

-^/^pWpWrfV. (11) 

The variation of the first term on the right-hand side can 
be written as 

5^ni£., = ^ [riiSsi + 5n^{£i + fc^)]. (12) 

i i 

In metallic systems the occupancy can change with an ad- 
ditional potential while in the non-metallic systems the 
occupation number remains unchanged for small pertur- 
bation. Thus, for the non- metallic systems one can write 



d'^rnEi = ^n^fci 



(13) 



Having used the perturbation series one obtains 
8E^ = ^ n.„fct, = ^{v\8U + We + ^p^c \ v) 

V 

^ {v\5U + 5Ve + Spxe \c){c\6U + We + ^Pxc \ v) 
+ Q{6U + 6Ve+6pxe). (14) 

where v and c denote the occupied and unoccupied bands, 
respectively. The Q stands for higher terms. The first 
term can be expressed as 

^(W I 5U + 5Ve + Spxe I v) = ^(v I 5U I v) 

V V 

+ f p''{r)SVeir)d^r+ f p"{r)Spxeir)d^r. (15) 



Having used the expression for the variation of the 
Coulombic potential 



K(r) = 



Sp{r') ,3 > 



r — r' 



(16) 



the difference of the second term of Eq. ( [ll|) and the sec- 
ond term of Eq. ( |l5| ) can be written as 

p^ir)SVeir) - lsipir)Veir)) ^ -\ / ^^^^^dV. 



r — r' 



(17) 



Having expressed the exchange and correlation potential 
from Eq.(pX|), the difference of the third term of Eq.dl^ 
and the third term of Eq. ( |l5| ) can be written as 

p°5pxe-5 [-^PP] =^ [p°^xc- {p- p°) p-^ 

(18) 



For the first order variation we can write 

and the second order variation can be written as 
6^ 



(19) 



5p2 



(A.,-(.-/).^)^ = -^. (20) 



where the subscript means that variations are taken for 
p — p° . The variation of E^i can be finaUy rewritten into 
the form 

V 

{v\6U + 6Vc + 6p^c \c){c\6U + 5V^ + 5^^^ \ v) 



E 

v,c 

1 f S p{r)6p(r') 3 3 , 
-; — a ra r 



C71 Of, 



2 ./ \r — r 
Se^dr) 



Spii^j / 



6p{r)dp{r)d^r 



e{SU + SVc + dp.,c)- 



(21) 



For the first order variation of the exchange-correlation 
potential p,xc we can write 



'^M.clr) =('%#') Spir). 



Spir) Jo 



(22) 



As it was already mentioned one needs to calculate the 
second derivatives of energy E^i 



d^E, 



r\ 11,0 r^ l/.h. 



(23) 



u=0 



and therefore from Eq. (Elf) is obvious that only the linear 
variation of the electron density is needed. 



IV. VARIATION OF ELECTRON DENSITY 

The charge density p{r) is calculated from 

p(r) = ^n,^*(r)V'.(r), (24) 

i 

where rii is the occupation number and tpi are one- 
electron wavef unctions, the variation of p(r) is expressed 
as 

Sp = ^(n, + Sn^M* + S^:M + H^) - n,^°*^l 

(25) 

Having supposed that in the non-metallic systems the 
bui ~ one can write 



bp^Y. ^^^'l'TH^ + HXi^i + brM^)- (26) 



The variation of the wave function 



(27) 



and the variation of the electron density can be written 



i jk 



(28) 



The variation of the siave function is obtained from the 
Lippmann-Schwingert3 equation 



where Go{e) is the Green function 



Go{e) = J2- 



e~ Si 



(29) 



(30) 



and 5V^^ is the variation of the total Kohn-Sham po- 
tential 

5V^'^ ^5U + 5Vc + Sfi^c- (31) 

With the use of Eqs.( [lq , p^ ) it can be written as 

SV^' = SU+fp^d^r'+(^J^) Sp{r). (32) 
J |r-r'| V Spi^) Jo 

The elements Uij are then expressed as 

a,, = {^p° I S^,) = (V'O I Go{e.,)SVKs \ i>°). (33) 



Having combined equations (E8f) and ( |33| ) and applying 
operator 9!J ^ which contains the limit to the zero per- 
turbation (Eq.(0)) we obtain integral equation 

5^,oP(r)=E"^^°*W€W 

X (V.°|Go(£.)5,lJ(7 + T4 + M.c)|V'°) 

= ^n,^r(r)^0(r)(VA" | Go{e.)dl^U \ V'") 



»j 



+ ^n,Vr(r)V°(r)(V°|Go(e.) 



a r + ' 



(34) 



d'pir) 



v.°) 



\r-r'\ V 6pir) J, ^■■' 

which can be rewritten into form 

5^,aP(r) = dJ^^Mr) (35) 

''A{v,v')r{r',r")d^r'dlMr")d'r", 



where a and A depend only on eigenfunctions and eigen- 
values of the unperturbed state. The quantities from 
previous equation can be expressed as 



a^,o^(r)-E' 



-e(r)V^°W(^,"|5^,,t/ |^°) 



«j 



Kv,v')^Y.- 



(36) 



:^^f(r)V'"(r)^;0(r')^°(r') (37) 



^,3 



r(r,r') 



5{t' 



5p{v") 



(38) 



After expansion into plane waves the integral equation 
( ^ ) yields the set of algebraic equations. By solving this 
set we obtain 9H ^p which are used for calculation of the 
dynamical matrix. 

The contribution of the variation of the Ef,i to the dy- 
namical matrix is given by the sum of contributions from 
four terms on the right-hand side of the Eq. (^), i.e. 



^a'a^lq,^ 



D^h (q,H- 



(39) 



The contribution of the fifth term on the right-hand side 
of the Eq. (|l|) is zero because 5«*„a;^_^e((5[/ + 6Vc + 

SlJixc) = 0. 

To calculate the term D)^^ (q, ^v) is quite easy because 
it doesn't depend on the variation of the charge density 
5p but on the variation of the pseudopotential 5U only. 
According to Eqs. (pT]) one can write 



(40) 



^o/(q>M^) = E(^'ki5ra^.%c^i«-k). 



More complicated is to calculate D'^u (q, /ijy) . According 



"^c.k+q 



to the Eq. (El^) one can write 






Dt'^^iq^H 



I ^v.k 
v,c,k 



^c.k+q 

x(«,k|a«:JC/-KK + A*:.c)|c,k + q) 
x(c,k-Hq|9^_^(f/ + K + M.c)|^;,k). 



(41) 



The sum over v and c in Eq.(M^) runs over occupied 



ei,3/ 



and unoccupied bands respectively. The D'^U (q, jmi/) and 



£)^i (q, ^Lv) can be written in the forms 



-»eZ,3 



^Ih (q' ^^^) 



!2>W«£f«,=,,3,, „,, 



r - r' 



and 



Dthi^.y-^) 



a/3 



5p{v) 



df:Mr)dl^pPir)d^r. 

I 

(43) 



V. CALCULATIONS 

The calculations were performed with all quantities, 
panded into plane waves. The set of special h-pointaL 
with (i i i) mesh was used for integration over the Bril- 
louin zone. It corresponds to 10 k-pomis for F phonon 
and 128 h-points for phonon with general wave-vector. 
The plane-wave set was selected for each k in the first 
Brillouin zone such that among the reciprocal-lattice vec- 
tors G only those are selected which have a kinetic energy 
Ekin = ^(k -f G)^ less then Ecutoff = l8Ry- The sums 
over unoccupied states went over energies up to cutoff 
energy Ecutoff = ^Ry, which corresponds to 350 unoccu- 
pied bands approximately. 



VI. RESULTS 

The method presented in this paper was applied to 
silicon crystal with lattice constant a = 0.543 nra. Both 
optical and acoustic branches of phonon frequencies in 
longitudinal and transverse modes were obtained for F, X 
and L phonons. Calculated frequencies together with 
values obtained by other authors and experimental values 
are summarized in Table 1. The values are expressed in 
THz. 

TABLE I. The calculated and experimental values of 
phonon frequencies in the silicon crystal. 





r 


L 


X 




(0,0,0) 


(- - -) 


(1,0,0) 


TO^alc'' 


15.9 


15.8 


15.1 


TO.alc'' 


15.6 


— 


— 


TOcalc" 


15.2 


— 


13.5 


TOcalc" 


15.7 


— 


— 


TOexpt' 


15.5 


15.0 


14.2 


LOcalc'^ 


15.9 


14.5 


14.6 


LOcalc'^ 


- 


- 


12.1 


LOexpt° 


15.5 


12.8 


12.5 


TA,,i,^ 


0.0 


6.2 


7.3 


TA^^iJ 


- 


3.0 


4.1 


TA^xpt" 


0.0 


3.6 


4.4 


LA.^ic'' 


0.0 


13.8 


14.6 


^^expt 


0.0 


11.7 


12.5 



^Phonons calculated in this work. 
^ Reference [Jl5 1 . 
^ Reference [ 



'^ Reference [ 
"Reference [ 
^ Reference [ 19 



VII. DISCUSSION 

Having used the presented method a quite good agree- 
ment of the calculated frequencies of optical phonons 
in silicon crystal and experimental ones was reached. 
We systematically increased the number of unoccupied 
states taken into account and we have found that the 
increase over ^ 300 bands almost doesn't change the re- 
sults. Thus, the contribution from states with energy 
£i > Scutoff doesn't appear essential. The biggest indi- 
vidual contributions to the dynamical matrix come from 
unoccupied states with energy near the Fermi energy. 
But the number of these states is small, the density of 
states is approximately g{e) ~ e^", and so it can be said 
that both dynamical matrix and phonon frequencies are 
integral quantities where the influence of one state is not 
essential. Any basis set for sums over unoccupied states is 
possible but the eigenfunctions of the Kohn-Sham hamil- 
tonian represent the natural choice . 

The presented method does not need time consuming 
calculations of inverse matrix as in the inuerse dielectric 
matrix method. Method of Baroni et ala is based on 
the self-consistent scheme for perturbed eigenfunctions 
while our method is based on the self-consistent equa- 
tion (B5[) for the variation of the electron density and 
only energies and eigenfunctions of unperturbed states 
are needed. King-Smith and Needa3 used the Hellmann- 
Feynman forces on all atoms in distorted crystal to con- 
strruct the dynamical matrix. The method of Gonze et 
alB is based on variations of the DFT total energy with 
respect to the first-order perturbations of the wave func- 
tions. In this paper, the second derivatives of energy ap- 
proach is used. The advantage of the presented method 
is that in contradiction to the both, Baroni et al. method 
and King-Smith and Needs method, the iteration up to 
the self-consistency is needed only once for the unper- 
tubed state. In contradiction to the method of Gonze 
et al., we don't need any additional minimization of the 
total energy with respect to the first-order variations of 
the wave function for each phonon, because all quantities 
are calculated from quantities of unperturbed state. 

Because of avoiding the iterations and minimization 
the presented method is very suitable for massive paral- 
lelization of the computer code. 



VIII. CONCLUSION 

The method for ab initio calculation of frequency and 
polarization vector of phonons with an arbitrary wave- 
vector was presented. The density functional theory 
within the pseudopotential framework was used. As in 
other linear response methods, the presented one uses the 
perturbation theory. In contradiction to the other meth- 
ods the presented one is based upon a direct calculation 
of the dynamical matrix via second derivatives of the 



total crystal energy and no calculation of the Hellman- 
Feynman forces is needed. 

The presented method does not need time consum- 
ing calculations of inverse matrix as in the inverse di- 
electric matrix method. Another advantage is that in 
contradistinction to thjs. both Baroni et al.a method and 
King-Smith and Necdaa method, no iteration up to self- 
consistency for individual phonon is needed. Similarly, 
we don't need any total energy minimization as in the 
case of method of Rignanese, Michenaud and GonzcQ. 
The computing time saving due to the parallelization of 
the computer code is significantly high for the presented 
method. 
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